module analyse_lattice
use my_tracking_code
implicit none
real(dp) :: lam =2.0d0,fix0(nvd) ,fixexp(nvd)
integer, parameter :: nsa=5
integer  nmap(nsa),mfexperiment(nsa)
real(dp) :: thick_lam=3.e-40,tune0
character*16 expfile(nsa)
character*16 expfile_ai(nsa)
character*16 expfile_ai_r(nsa)
integer, parameter:: nre=6
integer  limc
COMPLEX(DP) AM(NRE,0:2)
real(dp) tw(3)
logical :: cs=.false.

type my_data
 integer nturn,nsa,nt(nsa)
 real(dp) ray(nsa,2,0:2000)
 real(dp) rayn(nsa,2,0:2000)
 real(dp) r2(nsa,0:2000)
 real(dp) tune(nsa),n(nsa)
 complex(dp) fm(nre,nsa)
end type my_data

type(my_data) experiment

contains

subroutine store_mydata(lat,ray,pos,nord)
implicit none
type(my_lattice) lat
integer nturn,pos,i,j,nord,mf,RR
real(dp) ray(nvd),lamt,r(lnv)
type(real_8) rayt(nvd)
type(normalform) n




expfile(1)='1.txt'
expfile(2)='2.txt'
expfile(3)='3.txt'
expfile(4)='4.txt'
expfile(5)='5.txt'
expfile_ai(1)='1ai.txt'
expfile_ai(2)='2ai.txt'
expfile_ai(3)='3ai.txt'
expfile_ai(4)='4ai.txt'
expfile_ai(5)='5ai.txt'
expfile_ai_r(1)='1ai_r.txt'
expfile_ai_r(2)='2ai_r.txt'
expfile_ai_r(3)='3ai_r.txt'
expfile_ai_r(4)='4ai_r.txt'
expfile_ai_r(5)='5ai_r.txt'


do i=1,nsa
  mfexperiment(i)=newfile
 open(unit=mfexperiment(i),file=expfile(i))
enddo

nmap=0
experiment%nsa=nsa
!lam=((abs(ray(1))+thick_lam)/abs(ray(1)))**(1.d0/dfloat(nsa))

lamt=lam
call alloc(rayt)

do i=1,nvd
 rayt(i)=ray(i)
enddo

do j=1,experiment%nsa

 do i=0,experiment%nturn
  experiment%ray(j,1,i)=rayt(1)
  experiment%ray(j,2,i)=rayt(2)
  call track(rayt,lat,pos)

 enddo

  lamt=lam*lamt
  do i=1,2
   rayt(i)=lamt*ray(i)
  enddo

 
enddo
if(nord==1) nord=6
limc=1
fix0=0.d0
fix0(3)=rayt(3)
call find_orbit_with_fpp(lat,pos,fix0,1)
call init(nord,1,0,0)
call alloc(n)
!goto 1000

rayt(1)=fix0(1)+(1.d0.mono.1)
rayt(2)=fix0(2)+(1.d0.mono.2)
rayt(3)=fix0(3)
  call track(rayt,lat,pos)
n=rayt
tune0=n%tune(1)


    tw(1)=(n%a_t%v(1).sub.'10')**2+(n%a_t%v(1).sub.'01')**2
    tw(2)=(n%a_t%v(2).sub.'10')**2+(n%a_t%v(2).sub.'01')**2
    tw(3)=(n%a_t%v(2).sub.'10')*(n%a_t%v(1).sub.'10')+(n%a_t%v(2).sub.'01')*(n%a_t%v(1).sub.'01')
    tw(3)=-tw(3)
    write(6,*) " beta ",tw(1)
    write(6,*) " gamma ",tw(2)
    write(6,*) " alpha ",tw(3)


n%a_t=n%a_t**(-1)
n%a_t=n%a_t.SUB.1

WRITE(6,*) "DO YOU WANT TO READ A^-1 "
READ(5,*) RR
IF(RR==1) THEN
 MF=NEWFILE
  open(unit=MF,file=expfile_ai(1))
  CALL DAINPUT(n%a_t,MF)
  CALL DAPRINT(n%a_t,mfexperiment(1))
  NMAP(1)=NMAP(1)+1

  MF=closefile
elseif(rr==-1) then
  n%a_t=1
  CALL DAPRINT(n%a_t,mfexperiment(1))
  NMAP(1)=NMAP(1)+1
else
  CALL DAPRINT(n%a_t,mfexperiment(1))
  NMAP(1)=NMAP(1)+1

ENDIF
r=0.d0
do j=1,experiment%nsa

 do i=0,experiment%nturn
 r(1)= experiment%ray(j,1,i)-fix0(1)
 r(2)= experiment%ray(j,2,i)-fix0(2)
 r(3)=fix0(3)


 r=n%a_t*r
 experiment%rayn(j,1,i)=r(1)
 experiment%rayn(j,2,i)=r(2)
 experiment%r2(j,i)=r(1)**2+r(2)**2

 enddo
 enddo
1000 continue
call kill(n)
call kill(rayt)

RETURN

fixexp=0.d0


do i=0,experiment%nturn


 fixexp(1)= experiment%ray(1,1,i)+fixexp(1)
 fixexp(2)= experiment%ray(1,2,i)+fixexp(2)

enddo

fixexp=fixexp/(experiment%nturn+1)


write(6,*) " Approximate centre ",fixexp(1:2)

do i=0,experiment%nturn
do j=1,experiment%nsa

 experiment%rayn(j,1,i)= experiment%ray(j,1,i)-fixexp(1)
 experiment%rayn(j,2,i)= experiment%ray(j,2,i)-fixexp(2)
 experiment%r2(j,i)=experiment%rayn(j,1,i)**2+experiment%rayn(j,2,i)**2

enddo
enddo



  mf=newfile
 open(unit=mf,file="centre.txt")
write(mf,*) fixexp(1:2)
  mf=closefile


end subroutine store_mydata


subroutine store_fake(nord)
implicit none
integer nord,i,mf(nsa),mf0,limct,j
type(damap) ai(nsa)

limct=limc+1
if(nord<6) nord=3
!nmap=0

mf0=newfile
 open(unit=mf0,file="centre.txt")
do i=1,limc
  mf(i)=newfile
 open(unit=mf(i),file=expfile_ai(limc))
enddo


call init(nord,1,0,0)
do i=1,limc; call alloc(ai(i));call dainput(ai(i),mf(i)); enddo;


read(mf0,*) fixexp(1:2)  
  
do i=0,experiment%nturn
do j=limct,limct

 experiment%rayn(j,1,i)= experiment%ray(j,1,i)-fixexp(1)
 experiment%rayn(j,2,i)= experiment%ray(j,2,i)-fixexp(2)
! experiment%r2(j,i)=experiment%rayn(j,1,i)**2+experiment%rayn(j,2,i)**2
enddo
enddo
write(6,*) " limc = ",limc,limct
call print(ai(limc),6)

do i=1,limc
call Trans_data(ai(i),limct) 
enddo 
  






do i=1,limc
  mf(i)=closefile
enddo
mf0=closefile

do i=1,limc; call kill(ai(i)); enddo;


limc=limct
end subroutine store_fake



subroutine close_first_order
implicit none
integer i,mf,k
type(damap) Ai,ai_t,a_t
type(dragtfinn) df
type(taylorresonance) hr
type(taylor) inv
real(dp) beta,gamma,alpha


do i=1,limc
 rewind mfexperiment(i)
enddo

call alloc(ai,ai_t,a_t)
call alloc(df)
call alloc(hr)
call alloc(inv)

do k=1,limc


ai_t=1

do i=1,nmap(k)
 call dainput(ai,mfexperiment(k))
 ai_t=ai.o.ai_t
enddo
write(6,*) " fix ",fixexp(1:2)
call print(ai_t,6)


a_t=ai_t**(-1)


    beta=(a_t%v(1).sub.'10')**2+(a_t%v(1).sub.'01')**2
    gamma=(a_t%v(2).sub.'10')**2+(a_t%v(2).sub.'01')**2
    alpha=(a_t%v(2).sub.'10')*(a_t%v(1).sub.'10')+(a_t%v(2).sub.'01')*(a_t%v(1).sub.'01')
    alpha=-alpha
    write(6,*) " beta ",beta
    write(6,*) " gamma ",gamma
    write(6,*) " alpha ",alpha

if(cs) then
 ai_t%v(1)=(1.d0/sqrt(beta))*(1.d0.mono.1)
 ai_t%v(2)=(sqrt(beta))*(1.d0.mono.2) + (alpha/sqrt(beta))*(1.d0.mono.1)
endif


if(k==limc) then
 mf=newfile
 open(unit=mf,file=expfile_ai(k))
 call print(ai_t,mf)
 mf=closefile
endif
enddo  ! k


!do i=1,limc
!mfexperiment(i)=closefile
!enddo

ai_t=1
do i=1,limc
mf=newfile
 open(unit=mf,file=expfile_ai(i))
 call dainput(ai,mf)

 ai_t=ai.o.ai_t
mf=closefile

call print(ai,6)
enddo


DF=ai_t
inv=(1.d0.mono.'2')+(1.d0.mono.'02')
inv=texp(df%nonlinear,inv)
!inv=inv*ai_t
hr=inv

mf=newfile
open(unit=mf,file=expfile_ai_r(limc))

! call print(df%linear,mf)
 call print(hr%cos,mf)
 call print(hr%sin,mf)
 
 mf=closefile





call kill(df)
call kill(ai,ai_t,a_t)
call kill(hr)
call kill(inv)

end subroutine close_first_order


subroutine tune_mydata
implicit none
integer i,j
real(dp) a
 
 experiment%tune=0.d0

do j=limc,limc

 do i=0,experiment%nturn-1

 call phasead_double(experiment%rayn(j,1,i),experiment%rayn(j,2,i),experiment%rayn(j,1,i+1),experiment%rayn(j,2,i+1),a)

 experiment%tune(j)=experiment%tune(j)+a

 enddo

  experiment%tune(j)= experiment%tune(j)/experiment%nturn
enddo



write(6,*) " tunes "

do j=limc,limc
 write(6,*) j,experiment%tune(j)
enddo
write(6,*) " ***********************  "

end subroutine tune_mydata

subroutine phasead_double(xb,yb,xa,ya,a)
implicit none
real(dp) xa,ya,xb,yb
real(dp) a


a=-atan2(ya,xa)+atan2(yb,xb)

a=a/2.0d0/pi

if(a<0.0d0) a=a+1.0d0
if(a>1.0d0) a=a-1.0d0



end subroutine phasead_double

subroutine res_mydata0
use imsl
implicit none
integer i,j,m(nre),m1(nre),m2(nre),dm,dmp,N1,IPATH,nl,j11
complex(dp) e,q0,d(0:1,nsa),q(nre,nre)
complex(dp)  x(nre),b(nre), ma(nre,nre)
type(pbresonance) hr
type(pbfield) h
type(damap) map
type(complextaylor) ct
type(taylor)v(2)
real(dp) xk
real(dp) rx0,lam2
cs=.true.
N1=experiment%nturn+1
! Write(6,*) " using ",n1,"data points"

experiment%fm=0.d0
experiment%tune(1)=tune0
write(6,*) "tune0= ",tune0

m2(1)=2
m1(1)=0
m(1)=m2(1)-m1(1)

j11=2
m2(j11)=1
m1(j11)=1
m(j11)=m2(j11)-m1(j11)



write(6,*)" Compute Fourier Modes "


nl=2

do j=1,1
do dm=1,nl
 experiment%fm(dm,j)=0.d0
 e=exp(-i_*twopi* experiment%tune(j)*m(dm))

 do i=0,experiment%nturn
  experiment%fm(dm,j)=experiment%fm(dm,j)+experiment%r2(j,i)*e**i/( experiment%nturn + 1 )
 enddo
   
enddo

enddo
do j=1,1

q=0.d0
do dm=1,nl
b(dm)=experiment%fm(dm,j)
 do dmp=1,nl
 if(dm/=dmp) then
  q(dm,dmp)=exp(i_*twopi* experiment%tune(j)*(m(dm)-m(dmp)))
  q(dm,dmp)=(1.d0-q(dm,dmp)**n1)/(1.d0-q(dm,dmp))/n1 
 else
  q(dm,dmp)=1.d0
 endif
 enddo
enddo
ma=0.d0
do dm=1,nl
 ma(dm,nl)=ma(dm,nl)+q(nl,dm)
do dmp=1,nl-1
 ma(dm,dmp)= ma(dm,dmp)+2*i_*m(dmp)*(q(dmp,dm)-q(nl,dm))
enddo
enddo

x=0.d0

IPATH=1
CALL dLSACG (nl, ma, nre, B, IPATH, X)
do dm=1,nl
!write(6,*) x(1),abs(x(1))
x(dm)=x(dm)/(experiment%rayn(j,1,0)-i_*experiment%rayn(j,2,0))**m2(dm)/(experiment%rayn(j,1,0)+i_*experiment%rayn(j,2,0))**m1(dm)
experiment%fm(dm,j)=x(dm)
!if(mod(nmap,50)==0) write(6,*) m(dm),x(dm),abs(x(dm))
write(6,*)m(dm),abs(x(dm))    !x(dm),
enddo
enddo

!e=(experiment%fm(1,2)-lam**2*experiment%fm(1,1))/(1.d0-lam**2)







do j=1,1


call alloc(h); call alloc(map); call alloc(hr); call alloc(ct);call alloc(v)
!j=1

 
xk=-1.d0
ct=0.d0
v(1)=1.d0.mono.1
v(2)=1.d0.mono.2

do dm=1,nl-1

 ct=ct+xk*experiment%fm(dm,j)*v(1)**m1(dm)*v(2)**m2(dm)
 ct=ct+xk* CONJG(experiment%fm(dm,j)) *v(1)**m2(dm)*v(2)**m1(dm)

enddo

hr%cos%h=ct%r
hr%sin%h=ct%i

h=hr

map=1
map=texp(h,map)

call Trans_data(map,j)

nmap(j)=nmap(j)+1
call print(map,mfexperiment(j))
call kill(h); call kill(map); call kill(hr);call kill(ct);call kill(v)

enddo


write(6,*)" ############################################# "


end subroutine res_mydata0

subroutine res_mydata1
use imsl
implicit none
integer i,j,m(nre),m1(nre),m2(nre),dm,dmp,N1,IPATH,nl,j11
complex(dp) e,q0,d(0:1,nsa),q(nre,nre)
complex(dp)  x(nre),b(nre), ma(nre,nre)
type(pbresonance) hr
type(pbfield) h
type(damap) map
type(complextaylor) ct
type(taylor)v(2)
real(dp) xk
real(dp) rx0,lam2

N1=experiment%nturn+1
! Write(6,*) " using ",n1,"data points"

experiment%fm=0.d0

m2(1)=3
m1(1)=0
m(1)=m2(1)-m1(1)

m2(2)=2
m1(2)=1
m(2)=m2(2)-m1(2)

m2(3)=4
m1(3)=0
m(3)=m2(3)-m1(3)

m2(4)=3
m1(4)=1
m(4)=m2(4)-m1(4)

!m2(5)=5
!m1(5)=0
!m(5)=m2(5)-m1(5)


j11=5
m2(j11)=1
m1(j11)=1
m(j11)=m2(j11)-m1(j11)



write(6,*)" Compute Fourier Modes "


nl=5

do j=1,1
do dm=1,nl
 experiment%fm(dm,j)=0.d0
 e=exp(-i_*twopi* experiment%tune(j)*m(dm))

 do i=0,experiment%nturn
  experiment%fm(dm,j)=experiment%fm(dm,j)+experiment%r2(j,i)*e**i/( experiment%nturn + 1 )
 enddo
   
enddo

enddo
do j=1,1

q=0.d0
do dm=1,nl
b(dm)=experiment%fm(dm,j)
 do dmp=1,nl
 if(dm/=dmp) then
  q(dm,dmp)=exp(i_*twopi* experiment%tune(j)*(m(dm)-m(dmp)))
  q(dm,dmp)=(1.d0-q(dm,dmp)**n1)/(1.d0-q(dm,dmp))/n1 
 else
  q(dm,dmp)=1.d0
 endif
 enddo
enddo
ma=0.d0
do dm=1,nl
 ma(dm,nl)=ma(dm,nl)+q(nl,dm)
do dmp=1,nl-1
 ma(dm,dmp)= ma(dm,dmp)+2*i_*m(dmp)*(q(dmp,dm)-q(nl,dm))
enddo
enddo

x=0.d0

IPATH=1
CALL dLSACG (nl, ma, nre, B, IPATH, X)
do dm=1,nl
!write(6,*) x(1),abs(x(1))
x(dm)=x(dm)/(experiment%rayn(j,1,0)-i_*experiment%rayn(j,2,0))**m2(dm)/(experiment%rayn(j,1,0)+i_*experiment%rayn(j,2,0))**m1(dm)
experiment%fm(dm,j)=x(dm)
!if(mod(nmap,50)==0) write(6,*) m(dm),x(dm),abs(x(dm))
write(6,*)m(dm),abs(x(dm))    !x(dm),
enddo
enddo

!e=(experiment%fm(1,2)-lam**2*experiment%fm(1,1))/(1.d0-lam**2)







do j=1,1


call alloc(h); call alloc(map); call alloc(hr); call alloc(ct);call alloc(v)
!j=1

 
xk=-1.d0
ct=0.d0
v(1)=1.d0.mono.1
v(2)=1.d0.mono.2

do dm=1,nl-1

 ct=ct+xk*experiment%fm(dm,j)*v(1)**m1(dm)*v(2)**m2(dm)
 ct=ct+xk* CONJG(experiment%fm(dm,j)) *v(1)**m2(dm)*v(2)**m1(dm)

enddo

hr%cos%h=ct%r
hr%sin%h=ct%i

h=hr

map=1
map=texp(h,map)

call Trans_data(map,j)

nmap(j)=nmap(j)+1
call print(map,mfexperiment(j))
call kill(h); call kill(map); call kill(hr);call kill(ct);call kill(v)

enddo


write(6,*)" ############################################# "


end subroutine res_mydata1

subroutine res_mydata2
use imsl
implicit none
integer i,j,m(nre),m1(nre),m2(nre),dm,dmp,N1,IPATH,nl,j11
complex(dp) e,q0,d(0:1,nsa),q(nre,nre)
complex(dp)  x(nre),b(nre), ma(nre,nre)
type(pbresonance) hr
type(pbfield) h
type(damap) map
type(complextaylor) ct
type(taylor)v(2)
real(dp) xk
real(dp) rx0,lam2

N1=experiment%nturn+1
! Write(6,*) " using ",n1,"data points"

experiment%fm=0.d0

m2(1)=4
m1(1)=1
m(1)=m2(1)-m1(1)

m2(2)=3
m1(2)=2
m(2)=m2(2)-m1(2)

m2(3)=5
m1(3)=0
m(3)=m2(3)-m1(3)


j11=4
m2(j11)=1
m1(j11)=1
m(j11)=m2(j11)-m1(j11)



write(6,*)" Compute Fourier Modes "


nl=4

do j=1,1
do dm=1,nl
 experiment%fm(dm,j)=0.d0
 e=exp(-i_*twopi* experiment%tune(j)*m(dm))

 do i=0,experiment%nturn
  experiment%fm(dm,j)=experiment%fm(dm,j)+experiment%r2(j,i)*e**i/( experiment%nturn + 1 )
 enddo
   
enddo

enddo
do j=1,1

q=0.d0
do dm=1,nl
b(dm)=experiment%fm(dm,j)
 do dmp=1,nl
 if(dm/=dmp) then
  q(dm,dmp)=exp(i_*twopi* experiment%tune(j)*(m(dm)-m(dmp)))
  q(dm,dmp)=(1.d0-q(dm,dmp)**n1)/(1.d0-q(dm,dmp))/n1 
 else
  q(dm,dmp)=1.d0
 endif
 enddo
enddo
ma=0.d0
do dm=1,nl
 ma(dm,nl)=ma(dm,nl)+q(nl,dm)
do dmp=1,nl-1
 ma(dm,dmp)= ma(dm,dmp)+2*i_*m(dmp)*(q(dmp,dm)-q(nl,dm))
enddo
enddo

x=0.d0

IPATH=1
CALL dLSACG (nl, ma, nre, B, IPATH, X)
do dm=1,nl
!write(6,*) x(1),abs(x(1))
x(dm)=x(dm)/(experiment%rayn(j,1,0)-i_*experiment%rayn(j,2,0))**m2(dm)/(experiment%rayn(j,1,0)+i_*experiment%rayn(j,2,0))**m1(dm)
experiment%fm(dm,j)=x(dm)
!if(mod(nmap,50)==0) write(6,*) m(dm),x(dm),abs(x(dm))
write(6,*)m(dm),abs(x(dm))    !x(dm),
enddo
enddo

!e=(experiment%fm(1,2)-lam**2*experiment%fm(1,1))/(1.d0-lam**2)







do j=1,1


call alloc(h); call alloc(map); call alloc(hr); call alloc(ct);call alloc(v)
!j=1

 
xk=-1.d0
ct=0.d0
v(1)=1.d0.mono.1
v(2)=1.d0.mono.2

do dm=1,nl-1

 ct=ct+xk*experiment%fm(dm,j)*v(1)**m1(dm)*v(2)**m2(dm)
 ct=ct+xk* CONJG(experiment%fm(dm,j)) *v(1)**m2(dm)*v(2)**m1(dm)

enddo

hr%cos%h=ct%r
hr%sin%h=ct%i

h=hr

map=1
map=texp(h,map)

call Trans_data(map,j)

nmap(j)=nmap(j)+1
call print(map,mfexperiment(j))
call kill(h); call kill(map); call kill(hr);call kill(ct);call kill(v)

enddo


write(6,*)" ############################################# "


end subroutine res_mydata2


subroutine Trans_data(m,j)
implicit none
type(damap) m
real(dp) r(lnv)
integer i,j 

r=0.d0

!do j=1,experiment%nsa

 do i=0,experiment%nturn
 r(1)= experiment%rayn(j,1,i)   !-fix0(1)
 r(2)= experiment%rayn(j,2,i)   !-fix0(2)
 r(3)=fix0(3)

 r=m*r
 experiment%rayn(j,1,i)=r(1)
 experiment%rayn(j,2,i)=r(2)
 experiment%r2(j,i)=r(1)**2+r(2)**2
 enddo


 
!enddo



end subroutine Trans_data


subroutine find_orbit_with_fpp(lat,pos,fix0,nturn)
implicit none
type(my_lattice) lat
integer, optional :: nturn
integer pos,i,k,ntu,l
real(dp) fix(nvd),fix0(nvd),eps,normb,norma,f(2)
type(real_8) x(nvd)
type(damap) Id
type(gmap) g
logical doit
eps=1.d-4
call init(1,1,0,0) 
call alloc(x); call alloc(id); call alloc(g,2);
doit=.true.
normb=1.d10

ntu=1
if(present(nturn)) ntu=nturn


fix=fix0
k=0
do while(doit) 
k=k+1
id=1
do i=1,2
 x(i)=fix(i)+id%v(i)
enddo
 X(3)=FIX(3)
do l=1,ntu
 call track(x,lattice,pos)
enddo
f(1)=x(1)
f(2)=x(2)
!write(6,*) f
!pause 123
do i=1,2
g%v(i)=x(i)-fix(i)-morph(id%v(i))
enddo

g=g.oo.(-1)


do i=1,2
fix(i)=(g%v(i).sub.0)+fix(i)
enddo
norma=abs(g%v(1).sub.'0')+abs(g%v(2).sub.'0')

if(norma<eps) then
 if(norma>=normb) then
  doit=.false.
 else
  normb=norma
 endif
endif

if(k>10000) then
write(6,*) " fixed point not found in 10000 iterations "
  goto 100
 endif
enddo

fix0=fix

100 call kill(x);
call kill(id);call kill(g);
end subroutine find_orbit_with_fpp

subroutine make_twis_table_with_A(lat,delta,filename)
implicit none
type(my_LATTICE)LAT
type(real_8) x(3)
type(damap) a1,map
type(normalform) Normal
real(dp) delta,fix(3),s
integer i,mf
character*(*) filename

mf=20 
open(unit=mf,file=filename(1:len_trim(filename)))

fix=0.d0; fix(3)=delta;
call find_orbit_with_fpp(lat,1,fix)

call init(1,1,0,0) 
call alloc(x); call alloc(map,a1); call alloc(normal);


do i=1,2
 x(i)=fix(i)+(1.0d0.mono.i)
enddo
x(3)=delta

call track(x,lat,1)

do i=1,2
 map%v(i)=x(i)
enddo

normal=map

a1=normal%a_t

do i=1,2
 x(i)=fix(i)+a1%v(i)
enddo
x(3)=delta

   lat%magnet(1)%twiss(beta)=(x(1).sub.'10')**2+(x(1).sub.'01')**2
   lat%magnet(1)%twiss(alpha)=-(x(1).sub.'10')*(x(2).sub.'10')-(x(1).sub.'01')*(x(2).sub.'01')
   lat%magnet(1)%twiss(gamma)=(x(2).sub.'10')**2+(x(2).sub.'01')**2


do i=1,lat%n

   call track(x,lat,i,i+1)
   
   lat%magnet(i+1)%twiss(beta)=(x(1).sub.'10')**2+(x(1).sub.'01')**2
   lat%magnet(i+1)%twiss(alpha)=-(x(1).sub.'10')*(x(2).sub.'10')-(x(1).sub.'01')*(x(2).sub.'01')
   lat%magnet(i+1)%twiss(gamma)=(x(2).sub.'10')**2+(x(2).sub.'01')**2


enddo

s=0.d0
do i=1,lat%n
 write(mf,100) i,s,lat%magnet(i)%twiss(beta:gamma)
s=s+lat%magnet(i)%L
enddo
 write(6,*) "  In front of magnet 1 "
 write(6,100) 1,0.d0,lat%magnet(1)%twiss(beta:gamma)
 write(6,*) "  Tracked value after last magnet"
 write(6,100) i,s,lat%magnet(lat%n+1)%twiss(beta:gamma)

100  FORMAT( I4, 4(1x,F8.4))

call kill(x); call kill(map,a1); call kill(normal);

close(mf)
end subroutine make_twis_table_with_A


subroutine make_twis_table_with_invariant(lat,delta,filename)
implicit none
type(my_LATTICE)LAT
type(real_8) x(3)
type(damap) a1,map
type(normalform) Normal
real(dp) delta,fix(3),s
integer i,mf
character*(*) filename
type(taylor) Invariant,circle

mf=20 
open(unit=mf,file=filename(1:len_trim(filename)))

fix=0.d0; fix(3)=delta;
call find_orbit_with_fpp(lat,1,fix)

call init(2,1,0,0) 
call alloc(x); call alloc(map,a1); call alloc(normal);call alloc(Invariant,circle);

circle=(1.d0.mono.'2')+(1.d0.mono.'02')

do i=1,2
 x(i)=fix(i)+(1.0d0.mono.i)
enddo
x(3)=delta

call track(x,lat,1)

do i=1,2
 map%v(i)=x(i)
enddo

normal=map

a1=normal%a_t

do i=1,2
 x(i)=fix(i)+a1%v(i)
enddo
x(3)=delta


invariant=circle*(a1**(-1))

   lat%magnet(1)%twiss(beta)=invariant.sub.'02'
   lat%magnet(1)%twiss(alpha)=(invariant.sub.'11')/2.d0
   lat%magnet(1)%twiss(gamma)=invariant.sub.'20'


do i=1,lat%n

   call track(x,lat,i,i+1)
   
a1%v(1)=x(1)
a1%v(2)=x(2)
invariant=circle*(a1**(-1))

   lat%magnet(i+1)%twiss(beta)=invariant.sub.'02'
   lat%magnet(i+1)%twiss(alpha)=(invariant.sub.'11')/2.d0
   lat%magnet(i+1)%twiss(gamma)=invariant.sub.'20'


enddo

s=0.d0
do i=1,lat%n
 write(mf,100) i,s,lat%magnet(i)%twiss(beta:gamma)
s=s+lat%magnet(i)%L
enddo
 write(6,*) "  In front of magnet 1 "
 write(6,100) 1,0.d0,lat%magnet(1)%twiss(beta:gamma)
 write(6,*) "  Tracked value after last magnet"
 write(6,100) i,s,lat%magnet(lat%n+1)%twiss(beta:gamma)

100  FORMAT( I4, 4(1x,F8.4))

call kill(x); call kill(map,a1); call kill(normal);call alloc(Invariant,circle);

close(mf)
end subroutine make_twis_table_with_invariant


subroutine fit_tune(delta,tune)
implicit none
integer i
type(real_8) x(nvd)
type(normalform) normal
type(damap) map
type(gmap) g
real(dp) fix(nvd),delta,tune
logical :: more
real(dp) eps,bef,aft
more=.true.

eps=1.0d-6

bef=1.d10
aft=0.d0


do while(more) 

fix=0.d0
fix(3)=delta
call find_orbit_with_fpp(lattice,1,fix)

call init(2,1,1,0)
call alloc(x); call alloc(normal);call alloc(map);call alloc(g,1);

do i=1,lattice%n
if(lattice%magnet(i)%type==1) then
 call make_it_knob(lattice%magnet(i)%k(2),1)
!  lattice%magnet(i)%k(2)%kind=3
! lattice%magnet(i)%k(2)%i=3
endif
enddo

x(1)=fix(1)+(1.d0.mono.1)
x(2)=fix(2)+(1.d0.mono.2)
x(3)=fix(3)

call track(x,lattice,1)

map%v(1)=x(1)
map%v(2)=x(2)

normal=map


g%v(1)=(normal%dhdj%v(2)<='00')-tune


g=g.oo.(-1)

aft=(g%v(1).sub.'0')
do i=1,lattice%n
if(lattice%magnet(i)%type==1) lattice%magnet(i)%k(2)%r=lattice%magnet(i)%k(2)%r+aft
enddo

if(aft<eps) then
  
 if(aft>=bef) more=.false.
  bef=aft

endif

call kill(x);call kill(normal);call kill(map);call alloc(g,1);
do i=1,lattice%n
if(lattice%magnet(i)%type==1) lattice%magnet(i)%k(2)%kind=1
enddo
!write(6,*) " More "
!read(5,*) more
enddo



end subroutine fit_tune


end module analyse_lattice
